Using EASI Sentinel-1 RTC Gamma0 data¶

This notebook demonstrates how to load and use Sentinel-1 Radiometric Terrain Corrected (RTC) Gamma0 data generated in EASI.

The processing uses SNAP-10 with a graph processing tool (GPT) xml receipe for RTC Gamma0 and its variants.

For most uses we recommend the smoothed 20 m product (sentinel1_grd_gamma0_20m). We can process the 10 m products (sentinel1_grd_gamma0_10m, sentinel1_grd_gamma0_10m_unsmooth) on request. Please also ask if you wish to trial other combinations of the parameters.

RTC Gamma0 product variants¶

sentinel1_grd_gamma0_20m sentinel1_grd_gamma0_10m sentinel1_grd_gamma0_10m_unsmooth
DEM
copernicus_dem_30 Y Y Y
Scene to DEM extent multiplier 3.0 3.0 3.0
SNAP operator
Apply-Orbit-File Y Y Y
ThermalNoiseRemoval Y Y Y
Remove-GRD-Border-Noise Y Y Y
Calibration Y Y Y
SetNoDataValue Y Y Y
Terrain-Flattening Y Y Y
Speckle-Filter Y Y N
Multilook Y Y N
Terrain-Correction Y Y Y
Output
Projection WGS84, epsg:4326 WGS84, epsg:4326 WGS84, epsg:4326
Pixel resolution 20 m 10 m 10 m
Pixel alignmentPixelIsArea = top-left PixelIsArea PixelIsArea PixelIsArea

Units and conversions¶

The sentinel1_grd_gamma0_* data are in Intensity units. Intensity can be converted to dB and amplitude, and vice-versa, with the following. Practical Xarray examples are given below.

Intensity to/from dB:

       dB = 10*log10(intensity)
intensity = 10**(dB/10)

Intensity to/from Amplitude:

intensity = amplitude * amplitude
amplitude = sqrt(intensity)

In this notebook we have two functions for xarray datasets/arrays, using numpy.

def intensity_to_db(x):
    return 10*numpy.log10(x)

def db_to_intensity(db):
    return numpy.pow(10, db/10.0)

Reference: https://forum.step.esa.int/t/what-stage-of-processing-requires-the-linear-to-from-db-command

Set up¶

Import required packages and functions¶

In [1]:
# Basic plots
%matplotlib inline
# import matplotlib.pyplot as plt
# plt.rcParams['figure.figsize'] = [12, 8]

# Common imports and settings
import os, sys, re
from pathlib import Path
from IPython.display import Markdown
import pandas as pd
pd.set_option("display.max_rows", None)
import xarray as xr
import numpy as np

# Datacube
import datacube
from datacube.utils.rio import configure_s3_access
from datacube.utils import masking
import odc.geo.xr
from dea_tools.plotting import display_map

# EASI tools
import git
repo = git.Repo('.', search_parent_directories=True).working_tree_dir  # This gets the current repo directory. Alternatively replace with the easi-notebooks repo path in your home directory
if repo not in sys.path: sys.path.append(repo)
from easi_tools import EasiDefaults, xarray_object_size
from easi_tools.notebook_utils import mostcommon_crs, initialize_dask, localcluster_dashboard, heading

# Data tools
import hvplot.xarray
# import geoviews
import cartopy.crs as ccrs

# Dask
import dask
from dask.distributed import Client, LocalCluster

EASI environment¶

This is for convenience only.

It allows this notebook to be used in any EASI deployment defined as part of the easi-notebooks repository. Please substitute with your own values if adapting the notebook for your own work.

In [2]:
easi = EasiDefaults()

family = 'sentinel-1'
# product = this.product(family)
product = 'sentinel1_grd_gamma0_20m'
display(Markdown(f'Default {family} product for "{easi.name}": [{product}]({easi.explorer}/products/{product})'))
Successfully found configuration for deployment "asia"

Default sentinel-1 product for "asia": sentinel1_grd_gamma0_20m

Dask and ODC¶

In [3]:
# Dask local cluster
cluster = LocalCluster(n_workers=4)
client = Client(cluster)
server = f'https://hub.{easi.domain}'  # Or replace if not using EasiDefaults
user = os.environ.get('JUPYTERHUB_SERVICE_PREFIX')  # Current user
dask.config.set({"distributed.dashboard.link": f'{server}{user}' + "proxy/{port}/status"})  # port is evaluated by dask
display(client)

# Or use Dask Gateway - this may take a few minutes
# cluster, client = initialize_dask(use_gateway=True, workers=4)
# display(client)

# ODC
dc = datacube.Datacube()
configure_s3_access(aws_unsigned=False, requester_pays=True, client=client);

# List measurements for the product
dc.list_measurements().loc[[product]]

Client

Client-8b7f9b90-9d74-11ef-80e0-62370a50615c

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: https://hub.asia.easi-eo.solutions/user/pag064/proxy/8787/status

Cluster Info

LocalCluster

c244db88

Dashboard: https://hub.asia.easi-eo.solutions/user/pag064/proxy/8787/status Workers: 4
Total threads: 8 Total memory: 24.00 GiB
Status: running Using processes: True

Scheduler Info

Scheduler

Scheduler-a68d846f-e14d-4ab4-a355-6ed821bca6b3

Comm: tcp://127.0.0.1:42299 Workers: 4
Dashboard: https://hub.asia.easi-eo.solutions/user/pag064/proxy/8787/status Total threads: 8
Started: Just now Total memory: 24.00 GiB

Workers

Worker: 0

Comm: tcp://127.0.0.1:35327 Total threads: 2
Dashboard: https://hub.asia.easi-eo.solutions/user/pag064/proxy/40293/status Memory: 6.00 GiB
Nanny: tcp://127.0.0.1:35483
Local directory: /tmp/dask-scratch-space/worker-3pvwqoly

Worker: 1

Comm: tcp://127.0.0.1:41249 Total threads: 2
Dashboard: https://hub.asia.easi-eo.solutions/user/pag064/proxy/34257/status Memory: 6.00 GiB
Nanny: tcp://127.0.0.1:35201
Local directory: /tmp/dask-scratch-space/worker-es2gp9yu

Worker: 2

Comm: tcp://127.0.0.1:35695 Total threads: 2
Dashboard: https://hub.asia.easi-eo.solutions/user/pag064/proxy/34863/status Memory: 6.00 GiB
Nanny: tcp://127.0.0.1:46253
Local directory: /tmp/dask-scratch-space/worker-6of1m908

Worker: 3

Comm: tcp://127.0.0.1:36525 Total threads: 2
Dashboard: https://hub.asia.easi-eo.solutions/user/pag064/proxy/34637/status Memory: 6.00 GiB
Nanny: tcp://127.0.0.1:35369
Local directory: /tmp/dask-scratch-space/worker-8gfz33sj
Out[3]:
name dtype units nodata flags_definition aliases add_offset scale_factor
product measurement
sentinel1_grd_gamma0_20m vh vh float32 intensity NaN NaN [gamma0_vh] NaN NaN
vv vv float32 intensity NaN NaN [gamma0_vv] NaN NaN
angle angle float32 degrees NaN NaN [local_incidence_angle, localincidenceangle] NaN NaN

Choose an area of interest¶

In [4]:
# Set your own latitude / longitude

# EASI defaults
# latitude = easi.latitude
# longitude = easi.longitude
# time = easi.time

# Australia GWW
# latitude = (-33, -32.6)
# longitude = (120.5, 121)
# time = ('2020-01-01', '2020-01-31')

# Example: PNG
latitude = (-4.26, -3.75)
longitude = (144.03, 144.74)
time = ('2020-01-01', '2020-05-31')

# Bangladesh
# latitude = (21.5, 23.5)
# longitude = (89, 90.5)
# time = ('2024-05-01', '2024-06-10')

# Vietnam
# epsg:32648
# latitude = (9.1, 9.9)
# longitude = (105.6, 106.4)
# time = ('2024-01-01', '2024-09-10')

display_map(longitude, latitude)
Out[4]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Load data¶

In [5]:
data = dc.load(
    product = product, 
    latitude = latitude,
    longitude = longitude,
    time = time,
    dask_chunks = {'latitude':2048, 'longitude':2048},      # Dask chunk size
    group_by = 'solar_day',                    # Group by day method
)

display(xarray_object_size(data))
display(data)
'Dataset size: 4.55 GB'
<xarray.Dataset> Size: 5GB
Dimensions:      (time: 45, latitude: 2550, longitude: 3550)
Coordinates:
  * time         (time) datetime64[ns] 360B 2020-01-03T08:39:20 ... 2020-05-3...
  * latitude     (latitude) float64 20kB -3.75 -3.75 -3.751 ... -4.26 -4.26
  * longitude    (longitude) float64 28kB 144.0 144.0 144.0 ... 144.7 144.7
    spatial_ref  int32 4B 4326
Data variables:
    vh           (time, latitude, longitude) float32 2GB dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
    vv           (time, latitude, longitude) float32 2GB dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
    angle        (time, latitude, longitude) float32 2GB dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
Attributes:
    crs:           EPSG:4326
    grid_mapping:  spatial_ref
xarray.Dataset
    • time: 45
    • latitude: 2550
    • longitude: 3550
    • time
      (time)
      datetime64[ns]
      2020-01-03T08:39:20 ... 2020-05-...
      units :
      seconds since 1970-01-01 00:00:00
      array(['2020-01-03T08:39:20.000000000', '2020-01-04T20:05:09.500000000',
             '2020-01-08T08:47:31.000000000', '2020-01-11T19:56:57.500000000',
             '2020-01-15T08:39:20.000000000', '2020-01-16T20:05:09.500000000',
             '2020-01-20T08:47:31.000000000', '2020-01-23T19:56:56.500000000',
             '2020-01-27T08:39:19.000000000', '2020-01-28T20:05:08.500000000',
             '2020-02-01T08:47:30.000000000', '2020-02-04T19:56:56.500000000',
             '2020-02-08T08:39:19.000000000', '2020-02-09T20:05:08.500000000',
             '2020-02-13T08:47:30.000000000', '2020-02-16T19:56:56.500000000',
             '2020-02-20T08:39:19.000000000', '2020-02-21T20:05:08.500000000',
             '2020-02-25T08:47:30.000000000', '2020-02-28T19:56:55.500000000',
             '2020-03-03T08:39:19.000000000', '2020-03-04T20:05:08.500000000',
             '2020-03-08T08:47:30.000000000', '2020-03-11T19:56:56.500000000',
             '2020-03-15T08:39:19.000000000', '2020-03-16T20:05:08.500000000',
             '2020-03-20T08:47:30.000000000', '2020-03-23T19:56:56.500000000',
             '2020-03-27T08:39:19.000000000', '2020-03-28T20:05:08.500000000',
             '2020-04-08T08:39:19.000000000', '2020-04-09T20:05:08.500000000',
             '2020-04-13T08:47:31.000000000', '2020-04-16T19:56:56.500000000',
             '2020-04-20T08:39:20.000000000', '2020-04-21T20:05:09.500000000',
             '2020-04-25T08:47:31.000000000', '2020-04-28T19:56:57.500000000',
             '2020-05-02T08:39:20.000000000', '2020-05-03T20:05:09.500000000',
             '2020-05-19T08:47:32.500000000', '2020-05-22T19:56:58.500000000',
             '2020-05-26T08:39:22.000000000', '2020-05-27T20:05:11.500000000',
             '2020-05-31T08:47:33.000000000'], dtype='datetime64[ns]')
    • latitude
      (latitude)
      float64
      -3.75 -3.75 -3.751 ... -4.26 -4.26
      units :
      degrees_north
      resolution :
      -0.0002
      crs :
      EPSG:4326
      array([-3.7501, -3.7503, -3.7505, ..., -4.2595, -4.2597, -4.2599])
    • longitude
      (longitude)
      float64
      144.0 144.0 144.0 ... 144.7 144.7
      units :
      degrees_east
      resolution :
      0.0002
      crs :
      EPSG:4326
      array([144.0301, 144.0303, 144.0305, ..., 144.7395, 144.7397, 144.7399])
    • spatial_ref
      ()
      int32
      4326
      spatial_ref :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]
      grid_mapping_name :
      latitude_longitude
      array(4326, dtype=int32)
    • vh
      (time, latitude, longitude)
      float32
      dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
      units :
      intensity
      nodata :
      nan
      crs :
      EPSG:4326
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 1.52 GiB 16.00 MiB
      Shape (45, 2550, 3550) (1, 2048, 2048)
      Dask graph 180 chunks in 1 graph layer
      Data type float32 numpy.ndarray
      3550 2550 45
    • vv
      (time, latitude, longitude)
      float32
      dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
      units :
      intensity
      nodata :
      nan
      crs :
      EPSG:4326
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 1.52 GiB 16.00 MiB
      Shape (45, 2550, 3550) (1, 2048, 2048)
      Dask graph 180 chunks in 1 graph layer
      Data type float32 numpy.ndarray
      3550 2550 45
    • angle
      (time, latitude, longitude)
      float32
      dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
      units :
      degrees
      nodata :
      nan
      crs :
      EPSG:4326
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 1.52 GiB 16.00 MiB
      Shape (45, 2550, 3550) (1, 2048, 2048)
      Dask graph 180 chunks in 1 graph layer
      Data type float32 numpy.ndarray
      3550 2550 45
    • time
      PandasIndex
      PandasIndex(DatetimeIndex([       '2020-01-03 08:39:20', '2020-01-04 20:05:09.500000',
                            '2020-01-08 08:47:31', '2020-01-11 19:56:57.500000',
                            '2020-01-15 08:39:20', '2020-01-16 20:05:09.500000',
                            '2020-01-20 08:47:31', '2020-01-23 19:56:56.500000',
                            '2020-01-27 08:39:19', '2020-01-28 20:05:08.500000',
                            '2020-02-01 08:47:30', '2020-02-04 19:56:56.500000',
                            '2020-02-08 08:39:19', '2020-02-09 20:05:08.500000',
                            '2020-02-13 08:47:30', '2020-02-16 19:56:56.500000',
                            '2020-02-20 08:39:19', '2020-02-21 20:05:08.500000',
                            '2020-02-25 08:47:30', '2020-02-28 19:56:55.500000',
                            '2020-03-03 08:39:19', '2020-03-04 20:05:08.500000',
                            '2020-03-08 08:47:30', '2020-03-11 19:56:56.500000',
                            '2020-03-15 08:39:19', '2020-03-16 20:05:08.500000',
                            '2020-03-20 08:47:30', '2020-03-23 19:56:56.500000',
                            '2020-03-27 08:39:19', '2020-03-28 20:05:08.500000',
                            '2020-04-08 08:39:19', '2020-04-09 20:05:08.500000',
                            '2020-04-13 08:47:31', '2020-04-16 19:56:56.500000',
                            '2020-04-20 08:39:20', '2020-04-21 20:05:09.500000',
                            '2020-04-25 08:47:31', '2020-04-28 19:56:57.500000',
                            '2020-05-02 08:39:20', '2020-05-03 20:05:09.500000',
                     '2020-05-19 08:47:32.500000', '2020-05-22 19:56:58.500000',
                            '2020-05-26 08:39:22', '2020-05-27 20:05:11.500000',
                            '2020-05-31 08:47:33'],
                    dtype='datetime64[ns]', name='time', freq=None))
    • latitude
      PandasIndex
      PandasIndex(Float64Index([            -3.7501,             -3.7503,             -3.7505,
                                -3.7507,             -3.7509,             -3.7511,
                                -3.7513,             -3.7515,             -3.7517,
                                -3.7519,
                    ...
                     -4.258100000000001,             -4.2583,  -4.258500000000001,
                                -4.2587,  -4.258900000000001,             -4.2591,
                    -4.2593000000000005,             -4.2595, -4.2597000000000005,
                                -4.2599],
                   dtype='float64', name='latitude', length=2550))
    • longitude
      PandasIndex
      PandasIndex(Float64Index([          144.0301,           144.0303, 144.03050000000002,
                              144.0307,           144.0309,           144.0311,
                    144.03130000000002,           144.0315,           144.0317,
                              144.0319,
                    ...
                              144.7381,           144.7383, 144.73850000000002,
                              144.7387,           144.7389,           144.7391,
                    144.73930000000001,           144.7395,           144.7397,
                              144.7399],
                   dtype='float64', name='longitude', length=3550))
  • crs :
    EPSG:4326
    grid_mapping :
    spatial_ref

Conversion and helper functions¶

In [6]:
# These functions use numpy, which should be satisfactory for most notebooks.
# Calculations for larger or more complex arrays may require Xarray's "ufunc" capability.
# https://docs.xarray.dev/en/stable/examples/apply_ufunc_vectorize_1d.html
#
# Apply numpy.log10 to the DataArray
# log10_data = xr.apply_ufunc(np.log10, data)

def intensity_to_db(da: 'xarray.DataArray'):
    # TODO:
    # xx = da.where(da > 0, np.nan)  # Set values <= 0 to NaN
    return 10*np.log10(da)

def db_to_intensity(ds: 'xarray'):
    return np.pow(10, ds/10.0)


def make_image(ds: 'xarray', frame_height=300, **kwargs):
    """Return a Holoviews object that can be displayed or combined"""
    # TODO select spatial dim names from the given xarray
    defaults = dict(
        cmap="Greys_r",
        x = 'longitude', y = 'latitude',
        rasterize=True,
        geo=True,
        frame_height=frame_height,
    )
    defaults.update(**kwargs)
    return ds.hvplot.image(**defaults)


def select_valid_time_layers(ds: 'xarray', percent: float):
    """Select time layers that have at least a given percentage of valid data (e.g., >=5%)

    Example usage:
      selected = select_valid_time_layers(ds, 0.05)
      filtered == ds.sel(time=selected)
    """
    # TODO select spatial dim names from the given xarray
    return ds.count(dim=['latitude','longitude']).values / (data.sizes['latitude']*data.sizes['longitude']) >= percent
In [7]:
# Optional time layer filter

selected = select_valid_time_layers(data.vv, 0.05)
data = data.sel(time=selected).persist()
In [8]:
db_data = intensity_to_db(data).persist()  # TODO: Apply only to selected bands

Plot the data¶

In [9]:
# A single time layer for VV and VH, with linked axes

vvplot = make_image(data.vv.isel(time=0), clim=(0, 0.5), title=f'VV ({data.time.dt.strftime("%Y-%m-%d %H:%M:%S").values[0]})', clabel='Intensity')
vhplot = make_image(data.vh.isel(time=0), clim=(0, 0.1), title=f'VH ({data.time.dt.strftime("%Y-%m-%d %H:%M:%S").values[0]})', clabel='Intensity')
vvplot + vhplot
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
Out[9]:
In [10]:
# Make a dB plot

vvplot = make_image(db_data.vv.isel(time=0), clim=(-30, -3), title=f'VV ({data.time.dt.strftime("%Y-%m-%d %H:%M:%S").values[0]})', clabel='dB')
vhplot = make_image(db_data.vh.isel(time=0), clim=(-30, -1), title=f'VH ({data.time.dt.strftime("%Y-%m-%d %H:%M:%S").values[0]})', clabel='dB')
vvplot + vhplot
Out[10]:
In [11]:
# Subplots for each time layer for VV, with linked axes

make_image(db_data.vh, clim=(-30, -3), robust=True).layout().cols(4)
Out[11]:
In [12]:
db_data.vh.plot.hist(bins=np.arange(-30, 0, 1), facecolor='red')
Out[12]:
(array([3.1075960e+06, 4.0301090e+06, 4.4976790e+06, 3.9098620e+06,
        2.4765660e+06, 1.1831090e+06, 4.7941300e+05, 2.2565400e+05,
        1.4956500e+05, 1.3116700e+05, 1.4566300e+05, 2.0239700e+05,
        3.6227400e+05, 8.3570200e+05, 2.4834020e+06, 9.4824540e+06,
        2.7861116e+07, 3.0549379e+07, 8.3279060e+06, 1.0488690e+06,
        3.9038000e+05, 1.7312700e+05, 5.6280000e+04, 1.3022000e+04,
        2.4230000e+03, 4.2600000e+02, 1.0300000e+02, 5.9000000e+01,
        5.0000000e+01]),
 array([-30., -29., -28., -27., -26., -25., -24., -23., -22., -21., -20.,
        -19., -18., -17., -16., -15., -14., -13., -12., -11., -10.,  -9.,
         -8.,  -7.,  -6.,  -5.,  -4.,  -3.,  -2.,  -1.]),
 <BarContainer object of 29 artists>)
No description has been provided for this image

Make an RGB image¶

For an RGB visualization we use the ratio between VH and VV.

In [13]:
# Add the vh/vv band
db_data['vh_vv'] = db_data.vh / db_data.vv

# Scale the measurements by their median so they have a similar range for visualization
med = db_data / db_data.median(dim=['latitude','longitude'])
med.persist()
Out[13]:
<xarray.Dataset> Size: 2GB
Dimensions:      (time: 12, latitude: 2550, longitude: 3550)
Coordinates:
  * time         (time) datetime64[ns] 96B 2020-01-04T20:05:09.500000 ... 202...
  * latitude     (latitude) float64 20kB -3.75 -3.75 -3.751 ... -4.26 -4.26
  * longitude    (longitude) float64 28kB 144.0 144.0 144.0 ... 144.7 144.7
    spatial_ref  int32 4B 4326
Data variables:
    vh           (time, latitude, longitude) float32 435MB dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
    vv           (time, latitude, longitude) float32 435MB dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
    angle        (time, latitude, longitude) float32 435MB dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
    vh_vv        (time, latitude, longitude) float32 435MB dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
xarray.Dataset
    • time: 12
    • latitude: 2550
    • longitude: 3550
    • time
      (time)
      datetime64[ns]
      2020-01-04T20:05:09.500000 ... 2...
      units :
      seconds since 1970-01-01 00:00:00
      array(['2020-01-04T20:05:09.500000000', '2020-01-16T20:05:09.500000000',
             '2020-01-28T20:05:08.500000000', '2020-02-09T20:05:08.500000000',
             '2020-02-21T20:05:08.500000000', '2020-03-04T20:05:08.500000000',
             '2020-03-16T20:05:08.500000000', '2020-03-28T20:05:08.500000000',
             '2020-04-09T20:05:08.500000000', '2020-04-21T20:05:09.500000000',
             '2020-05-03T20:05:09.500000000', '2020-05-27T20:05:11.500000000'],
            dtype='datetime64[ns]')
    • latitude
      (latitude)
      float64
      -3.75 -3.75 -3.751 ... -4.26 -4.26
      units :
      degrees_north
      resolution :
      -0.0002
      crs :
      EPSG:4326
      array([-3.7501, -3.7503, -3.7505, ..., -4.2595, -4.2597, -4.2599])
    • longitude
      (longitude)
      float64
      144.0 144.0 144.0 ... 144.7 144.7
      units :
      degrees_east
      resolution :
      0.0002
      crs :
      EPSG:4326
      array([144.0301, 144.0303, 144.0305, ..., 144.7395, 144.7397, 144.7399])
    • spatial_ref
      ()
      int32
      4326
      spatial_ref :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]
      grid_mapping_name :
      latitude_longitude
      array(4326, dtype=int32)
    • vh
      (time, latitude, longitude)
      float32
      dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
      Array Chunk
      Bytes 414.39 MiB 16.00 MiB
      Shape (12, 2550, 3550) (1, 2048, 2048)
      Dask graph 48 chunks in 1 graph layer
      Data type float32 numpy.ndarray
      3550 2550 12
    • vv
      (time, latitude, longitude)
      float32
      dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
      Array Chunk
      Bytes 414.39 MiB 16.00 MiB
      Shape (12, 2550, 3550) (1, 2048, 2048)
      Dask graph 48 chunks in 1 graph layer
      Data type float32 numpy.ndarray
      3550 2550 12
    • angle
      (time, latitude, longitude)
      float32
      dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
      Array Chunk
      Bytes 414.39 MiB 16.00 MiB
      Shape (12, 2550, 3550) (1, 2048, 2048)
      Dask graph 48 chunks in 1 graph layer
      Data type float32 numpy.ndarray
      3550 2550 12
    • vh_vv
      (time, latitude, longitude)
      float32
      dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
      Array Chunk
      Bytes 414.39 MiB 16.00 MiB
      Shape (12, 2550, 3550) (1, 2048, 2048)
      Dask graph 48 chunks in 1 graph layer
      Data type float32 numpy.ndarray
      3550 2550 12
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2020-01-04 20:05:09.500000', '2020-01-16 20:05:09.500000',
                     '2020-01-28 20:05:08.500000', '2020-02-09 20:05:08.500000',
                     '2020-02-21 20:05:08.500000', '2020-03-04 20:05:08.500000',
                     '2020-03-16 20:05:08.500000', '2020-03-28 20:05:08.500000',
                     '2020-04-09 20:05:08.500000', '2020-04-21 20:05:09.500000',
                     '2020-05-03 20:05:09.500000', '2020-05-27 20:05:11.500000'],
                    dtype='datetime64[ns]', name='time', freq=None))
    • latitude
      PandasIndex
      PandasIndex(Float64Index([            -3.7501,             -3.7503,             -3.7505,
                                -3.7507,             -3.7509,             -3.7511,
                                -3.7513,             -3.7515,             -3.7517,
                                -3.7519,
                    ...
                     -4.258100000000001,             -4.2583,  -4.258500000000001,
                                -4.2587,  -4.258900000000001,             -4.2591,
                    -4.2593000000000005,             -4.2595, -4.2597000000000005,
                                -4.2599],
                   dtype='float64', name='latitude', length=2550))
    • longitude
      PandasIndex
      PandasIndex(Float64Index([          144.0301,           144.0303, 144.03050000000002,
                              144.0307,           144.0309,           144.0311,
                    144.03130000000002,           144.0315,           144.0317,
                              144.0319,
                    ...
                              144.7381,           144.7383, 144.73850000000002,
                              144.7387,           144.7389,           144.7391,
                    144.73930000000001,           144.7395,           144.7397,
                              144.7399],
                   dtype='float64', name='longitude', length=3550))
In [14]:
# RGB plot using a DEA Tools function
# https://github.com/GeoscienceAustralia/dea-notebooks
# https://github.com/digitalearthafrica/deafrica-sandbox-notebooks/

sys.path.append('/home/jovyan/dea-notebooks/Tools')
from dea_tools.plotting import rgb

rgb(med, bands=['vh','vv','vh_vv'], col="time")
/env/lib/python3.10/site-packages/matplotlib/cm.py:494: RuntimeWarning: invalid value encountered in cast
  xx = (xx * 255).astype(np.uint8)
No description has been provided for this image
In [ ]:
# Experimental - Use Holoviews

# Create an RGB array, and persist it on the dask cluster
rgb_ds = xr.concat([med.vv, med.vh, med.vh_vv], 'channel').rename('rgb').to_dataset().persist()

# Plot the RGB
rgb_plot = rgb_ds.hvplot.rgb(
    bands='channel',
    groupby='time', rasterize=True,
    geo=True,
    title='RGB', frame_height=500,
)

rgb_plot  # + vv_plot + vh_plot

Export to Geotiffs¶

Recall that to write a dask dataset to a file requires the dataset to be .compute()ed. This may result in a large memory increase on your JupyterLab node if the area of interest is large enough, which in turn may kill the kernel. If so then skip this step, choose a smaller area or find a different way to export data.

In [ ]:
# Make a directory to save outputs to
target = Path.home() / 'output'
if not target.exists(): target.mkdir()

def write_band(ds, varname):
    """Write the variable name of the xarray dataset to a Geotiff files for each time layer"""
    for i in range(len(ds.time)):
        date = ds[varname].isel(time=i).time.dt.strftime('%Y%m%d').data
        single = ds[varname].isel(time=i).compute()
        write_cog(geo_im=single, fname=f'{target}/example_sentinel-1_{varname}_{date}.tif', overwrite=True)
        
write_band(data, 'vv')
write_band(data, 'vh')
# write_band(rgb_da, 'rgb')
In [ ]: